Multistable attractors in a network of phase oscillators with three-body interaction 
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Three-body interactions have been found in physics, biology, and sociology. To investigate their 
effect on dynamical systems, as a first step, we study numerically and theoretically a system of phase 
oscillators with three-body interaction. As a result, an infinite number of multistable synchronized 
states appear above a critical coupling strength, while a stable incoherent state always exists for 
any coupling strength. Owing to the infinite mult ist ability, the degree of synchrony in asymptotic 
state can vary continuously within some range depending on the initial phase pattern. 
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Interaction among particles or elements in classical 
mechanics, electromagnetism, and many other fields of 
physics is often modeled by two-body interaction. De- 
scription by the linear superposition of two-body inter- 
actions has allowed us to predict the future orbits of the 
planets and to design drug molecules that tightly bind 
to the target protein. However, it has been revealed that 
the net interaction experienced by an element cannot be 
written as the linear superposition of the two-body inter- 
action in several systems, including physical systems PQ, 
social and economic systems [2], and neuronal networks 
[3]. A typical example is signal transmission from one 
neuron to another. The signals are mediated by the re- 
lease of neurotransmitters from synapses, and some neu- 
rotransmitters modulate the response of neurons to in- 
puts from other neurons (heterosynaptic plasticity) [4]. 
This modulation can be regarded as three-body inter- 
action, although synaptic transmission is conventionally 
modeled as a two-body interaction. To show what oc- 
curs in such neuronal networks with three-body interac- 
tions, we present a numerical simulation of a network 
of Hodgkin-Huxley neurons [5] with short-term heterosy- 
naptic plasticity (see [6 for details). In this model, the 
input from neuron j to i is modulated by the relative 
spike timing of neuron i and other neurons in this model. 
Figure [TJa) shows that this neuronal network exhibits 
multistability, in which the numbers of synchronized neu- 
rons at the steady state vary depending on the initial 
conditions [Fig. [ljb)]. This seems to be a novel behavior 
not observed in systems with only two-body interactions. 
However, this system is too complicated to show analyt- 
ically why this multistability arises. 

To analyze the neuronal networks with three-body in- 
teraction, we exploit the fact that neurons exhibit peri- 
odic firings in many cases. Periodic activities are ubiq- 
uitous in not only neuronal networks, but also phenom- 
ena studied in other fields of biology, including gene ex- 
pression in E. col% synchronous flashing of fireflies, and 




FIG. 1. (color online). Two examples of multistability aris- 
ing from three-body interaction, (a) Raster plot of the spikes 
in a network of the Hodgkin-Huxley neurons with short-term 
heterosynaptic plasticity (N = 500). Black, red, and blue 
dots represent the firings in the networks starting from differ- 
ent initial conditions. Neurons are sorted in ascending order 
of the firing rate. Firing rate of each neuron is shown in (b). 
(c) Time evolution of the order parameter R for three dif- 
ferent initial conditions in the phase-oscillator systems with 
non- uniform random coupling (N = 500). Phase distribu- 
tions of oscillators at t = and t = 1 , 000 are shown on the 
left and right, respectively. These two examples demonstrate 
that the same system can show different degrees of synchrony 
depending on the initial conditions. 



pedestrians' gait [7 . The behavior of these periodic ac- 
tivities is described by a form of phase oscillators in a 
quite general context [8]. However, three-body interac- 
tion among phase oscillators has not been studied yet. 
Since phase oscillators are simple enough to be analyt- 
ically tractable and structurally stable, theory of phase 
oscillators is a powerful tool in interpreting and elucidat- 
ing complicated experimental results in which three-body 
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interactions play an essential role. In this Letter, we thus 
examine the effect of three-body interaction on the dy- 
namics of globally-coupled phase oscillators. 

As a natural extension of the system of limit-cycle 
oscillators with two-body interaction, iV-oscillator sys- 
tem with two- and three-body interaction is described 
by X, = F,(X,) + E^V^X^X^X,), where Fi 

describes the dynamics of uncoupled oscillator i and 
Vijk is the phase coupling function. Two-body inter- 
action V^-(X^Xj) is then included as a special case of 
the three-body interaction VV^X^, Xj, X&). Using the 
phase reduction technique, we can describe the dynamics 
of oscillator i with one variable, phase fa. Thus, the dy- 
namics of the system of phase oscillators with three-body 
interaction is generally given by 

fa = Ui + y^r ijfc (^,fc), (i) 

where UJi is the natural frequency of oscillator z, fai = 
<pj — fa, and Tijk is the coupling function. 

We present one example in which typical novel features 
arising from three-body interactions can be seen: 

fa = + y^Jaij sin(0 ji + auj) + sm(2fai + a 2i j)] 

3 

+ °ijk sin(0^ + otsijk) cos(fai + a Aijk ), 

where ~ A/"(0.3, 0.01), c^fc ~ A/"(6,4), and 

OLiij,a 2 ij,ot3ijk,cx4ijk ~ A/"(0,0.09). Here j\f(fi,cr 2 ) de- 
notes the normal distribution with mean /i and variance 
a 2 . In all simulations throughout this Letter, the nat- 
ural frequencies are drawn from A/"(0,1). Some typical 
time evolutions of the order parameter R representing 
the degree of synchrony is shown in Fig. [ljc), in which 
the above system starts from different initial conditions. 
The order parameter R is defined by 

Rexp(iip) = ^exp(i^), (2) 
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where ip is the average phase associated with the or- 
der parameter. As illustrated in Fig. [ljc), the sys- 
tem starting from a completely uniform initial distribu- 
tion remains desynchronized, while the system with non- 
uniform initial distribution can go to various synchro- 
nized states in a similar way as in Fig. [TJa). Two numer- 
ical simulations shown in Fig. [I] suggest that the system 
containing three-body interaction can exhibit multistable 
behaviors in a structurally stable manner. 

To investigate these behaviors analytically, we here im- 
pose three assumptions which do not spoil the essence of 
the above dynamical behaviors. First we assume that the 
phase coupling functions are identical for all oscillators, 
that is, Tijktfaiifai) = Y^ji^faij/N 2 . Second, without 
loss of generality, we can assume that the phase coupling 



function is symmetric, that is, T((pji,fai) = T(fai,(j)ji), 
because replacing the asymmetric coupling T asym (x,y) 
with symmetric coupling T 8ym (x,y) = [T asyrri (x,y) + 
r a sym(^ 7 #)]/2 does not change the dynamics. The 
last assumption is that inverting the phases of oscil- 
lators inverts the sign of forces among them, that is, 
r(0ji,^>/ci) = —T(faj,fak). Although this seems a rather 
tight constraint, this antisymmetricity is a property 
of the classical two-body coupling function T(4>ji) — 
sin^. We confirmed that the system under these con- 
straints could exhibit the qualitatively same behavior as 
in Fig. [TJc). Finally, we note that, owing to the 2tt- 
periodicity, T can be approximated by the finite Fourier 
series r(^,^ fci ) = K 2 (smfa i +sm(pki)/2+K 2 '(sm2(p j i + 
sin20fci)/2 + Ks{sm<pji cos fai + cos <pji sin^), where 
K 2l K' 2 and K3 are constants. Thus, the dynamics of 
globally-coupled phase oscillators with this type of three- 
body coupling is given by 

fa = UJi + y^(if 2 sin fai + K' 2 sin 2(j) ji ) 
3 

+ Sin fa* COS ' ( 3 ) 

3,k 

We further simplify this model equation to make it an- 
alytically tractable. Using order parameter R and setting 
K 2 = 0, K 3 = 0, and K 3 = K, we obtain the equations 
of dynamics with pure three-body interaction, 

0i =uJi- KR 2 sin20i, (4) 

where 0i = fa — ip is the relative phase of oscillator i 
to the average phase ip [Eq. Q]. Because we are using 
a co-rotating frame, we may here assume that the aver- 
age phase ip is constant. We assume that the frequency 
of the average phase ip equals the mean of the distribu- 
tion g(uj) of the natural frequency, the standard normal 
distribution. This assumption simplifies the equations 
to be derived, and, in addition, the solution of the de- 
rived self-consistent equation fits substantially well with 
the numerical results, although this assumption may not 
hold in some cases. 

Numerical simulations of Eq. Q with N = 10, 000 
oscillators and K = 3 from three initial conditions are 
shown in Fig. [2|a). Order parameter R takes various 
values depending on the initial conditions. Synchronized 
and desynchronized states coexist in the same parameter 
region. The relationship between the natural frequency 
UJi and the phase fa is also shown in Fig. [2|b,c,d). Fig- 
ure |2|c) indicates that oscillators can be phase locked to 
two specific phases. Indeed, an oscillator with natural fre- 
quency Ui can be phase locked to 6i = \ arcsin , tt + 
\ arcsin jf^z, if —KR 2 < UJi < KR 2 . On the other hand, 
Fig.[2|d) shows that the system with the same parameter 
values can exhibit a completely desynchronized state. 

If all of the phase-locked oscillators are locked to 0i = 
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FIG. 2. (a) Time evolution of the order parameter R from 
three different initial conditions in the mean field model with 
N — 10, 000 and if — 3. (b,c,d) uii-fa relationship for differ- 
ent initial conditions at t = 10, 000. 



\ arcsin , R is given by 
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g(u) du 



nir/4 

= 2KR 2 / cos cos 20g(KR 2 sin 20) d<9 

J-tt/4 

= S(R,K), (5) 

where we used du/ d6 = 2KR 2 cos2t9, and assumed that 
the non-phase-locked oscillators do not contribute to the 
value of the order parameter because the distribution 
g(uj) of natural frequency is the standard normal dis- 
tribution. The self-consistent equation R = S(R, if) has 
a solution R = for any if. In addition, S"(0,if) = 
^-\ R=0 = suggests that this solution is stable. For 
some if, the self-consistent equation has a non-zero so- 
lution R = r > or two non-zero solutions > r± > 
[Fig. |3|[a)]. Equation g gives the order parameter of 
the system in which all of the phase-locked oscillators 
take Oi = ^arcsin jf^2- Oscillator z, however, can also 
be phase locked to Oi = 7r + \ arcsin . Defining n(0) 
as the number of oscillators phase locked to 0, we charac- 
terize the distribution of the phase-locked oscillators with 
the function q(0) = [n(0) - n(0 + 7r)]/[n(0) + n(0 + tt)]. 
Note that \q{&)\ < 1. Then, the order parameter of the 
system is given by 



R = 2KR 2 



tt/4 



tt/4 



q(0)C{0, R, if) dO = S[R, if, q(0)], 



(6) 

where C(0,R,K) = cos cos 20 g (KR 2 sin 2(9). 

The largest attainable R for coupling strength if is 
given by the largest solution T2 of Eq. (J5|, while the small- 
est attainable non-zero R for the coupling strength if is 
given by the minimum of the largest positive solution of 
the self-consistent equation Eq. ([6| over all possible re- 
alizations of q(0). If R = S[R, K,q(0)] has two non-zero 



solutions V2 > ri, there exists < a < 1 with which 
the largest non-zero solution of R = S[R,K,aq(0)] is 
smaller than r2 because S[R, if, aq(0)] = aS[R,K,q(0)] 
[Fig. |3ja)]. Hence, to obtain the lowest attainable R, we 
have to find q(0) with which R = S[R,K,q(0)\ has only 
one non-zero solution. In other words, we find the small- 
est r satisfying S[r,K,q(0)] = r and 5"[r, if, g(0)] = 1 in 
varying q(0), where 

d 



S'[r,K,q(0)] 



OR 



S[R,K,q(0)] 



J -tt/4 



4K> / q(0)W(0,r,K)C(0,r,K)dO, 



and 



W{0,r,K) 



T , 2 • n J{Kr 2 ^i20) 
g(Kr z sin 20) 

[Fig. (3^c)]. To this end, first we fix R to r and examine 
whether there exists a solution — 1 < q(0) < 1 of the 
equations S[r, if , g(0)] = r and 5"[r, if, g(0)] = 1. If it 
exists, there is a solution — 1 < qi(0) < 1 of equations 
S[r,if,g(0)] = r and S'[r,K,q(0)] = s ± < 1 [Fig. §b), 
blue line], and there is a solution — 1 < q2(0) < 1 of 
equations S[r,K,q(0)] = r and S'[r,K,q(0)] = s 2 > 1 
[Fig. |3jb), green line]. Conversely, if qi(0) and #2(0) 
are given, -1 < q(0) = uq 2 (0) + (1 - u)gi(0) < 1, 
where < u = < 1, is a solution of the equa- 

tions S[r, if, q(0)} = r and S"[r, if, g(0)] = 1 [Fig. ^p, red 
line]. Thus, the existence of qi(0) and #2(0) which satisfy 
S[r,K, qi (9)] = r, S'[r,K, qi (0)] < 1, S[r,K,q 2 (9)] = r, 
and S'[r,K,q2(0)] > 1 is a necessary and sufficient con- 
dition of the existence of the solution of S[r, if, q(0)] = r 
and S f [r, if, q{6)\ = 1. It is sufficient for us to calculate 
the maximum and the minimum of S'[r,K,q(Q)\ under 
the constraints S[r,K,q(Q)\ = r and \q(0)\ < 1. 

S[r, if , q(0)] and S'[r,K,q(Q)\ have the same domain 
of integration, and their integrands differ by a factor of 
W(t9,r, if). Hence, the maximum of S'[r,K,q(Q)\ under 
the conditions S[r, if , q(0)] — r and \q(0)\ < 1 is given 
by S"[r, if, q 2 (0)} where q 2 {0) = 2G[W(0, r, if) - w 2 ] - 1. 
Here 9(x) is the Heaviside function, and is set to 
satisfy S[r, if, ^2(0)] = ^- I n this case, the phase-locked 
oscillators are distributed according to n(0)/[n(0)+n(0 + 
7r)] = 0[W(0,r, if) — ^2]. In other words, first we ad- 
just W2 to set S[r, if, ^2(0)] = r [Fig. [3jd)] , and next we 
check whether S"[r, if , #2(0)] is larger than 1 [Fig. j3^e)]. 
In the same way, we vary w\ to set S[r, if, ^i(0)] = r, 
where 91 (0) = 26[^i — W(0, r, if )] — 1, and check whether 
S'[r,K,qi(9)] is smaller than 1. 

Thus, we have theoretically obtained the region of the 
order parameter which can be achieved by choosing suit- 
able initial conditions [Fig. [3jf), red line]. In this figure, 
the dots represent the data from numerical simulations 
(N = 10, 000). The theoretical results agree with the nu- 
merical ones, though several points with if < 3 lie out- 
side of the theoretically derived region. This discrepancy 
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may be because the system size is too small or because 
we assumed that the frequency of the average phase co- 
incides with the mean of the distribution g{uo). 

Finally, we should remark that interactions in real- 
world systems generally contain not only three-body but 
also two-body interactions. We thus examine the behav- 
ior of the system described by 



= UJj 



N 



^2> 



s'm (t)ji cos (f)ki- 



As Fig.[3jg) shows, as K% increases, the system first starts 
out exhibiting either a single synchronized or desynchro- 
nized state depending on K2. Then briefly, a small win- 
dow in which these two states are bistable, appears. Fi- 
nally, multistable synchronized states, or for smaller if 2, 
a coexistence of desynchronized and multistable synchro- 
nized states [Fig. [3^g), orange region] corresponding to 
the multistability shown in Fig. [2j appears. This implies 
that our theoretical result derived with pure three-body 
interaction is structurally stable and generic. 

In this Letter, we have examined the behavior of phase 
oscillators with three-body interactions. We have found 
that this system can take an infinite number of synchro- 
nized states in a structurally stable manner [Fig. (3^g)]. 
We have derived the range of the order parameter R that 
can be attained by varying the initial condition. Our re- 
sults are different from the chimera state [9], because in 
our model we can continuously control the order param- 
eter of the steady state by choosing the initial condition. 
In addition, our model system can be completely inco- 
herent even in the K — > 00 limit (cf. [E]). There remain 
several questions to be answered. Three-body interac- 
tions in real-world systems and their behavior should be 
compared to those of the present model. Neurophysio- 
logical experiments [11] have shown that some prefrontal 
neurons keep their level of activity for several seconds. 
It is believed that this persistent activity serves as work- 
ing memory by encoding an analog quantity in the firing 
rate of multistable neuronal networks. Our results sug- 
gest the possibility that working memory uses the degree 
of synchrony among neurons to encode an analog quan- 
tity. Finally, we should systematically investigate various 
types of coupling function and the dynamical behavior on 
complex networks [T2] . 
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